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Multiscale problems are computationally costly to solve by direct sim- 
ulation because the smallest scales must be represented over a domain 
determined by the largest scales of the problem. We have developed and 
analyzed new numerical methods for multiscale wave propagation follow- 
ing the framework of the heterogeneous multiscale method. The numerical 
methods couple simulations on macro- and microscales for problems with 
rapidly fluctuating material coefficients. The computational complexity of 
the new method is significantly lower than that of traditional techniques. 
We focus on HMM approximation applied to long time integration of one- 
dimensional wave propagation problems in both periodic and non-periodic 
medium and show that the dispersive effect that appear after long time 
is fully captured. 

1 Introduction 

There are many application of the wave propagation in heterogeneous media. 
Acoustics, elastic or electromagnetic wave propagation in composite material 
are typical examples. In seismic exploration an acoustic wave equation must 
be solved for highly fluctuating wave velocity over long time with large distur- 
bances. We consider a representable model problem: The scalar wave equation 
in a bounded domain with initial and boundary conditions, 



on a smooth domain C R d where A £ (x) is a symmetric, uniformly positive 
definite matrix. For simplicity we assume that u e satisfies periodic boundary 
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conditions. We assume that A £ has oscillations on a scale proportional to e <C 1. 
The solution of ([I]) will then be a sum of two parts: one coarse scale (macroscale) 
part, which is independent of e, and an oscillatory (microscale) part which is 
highly oscillating in both time and spatial directions on the scale e. These 
kinds of multiscale problems are typically very computationally costly to solve 
by traditional numerical techniques. The smallest scale must be well represented 
over a domain which is determined by the largest scales. However most often 
one is only interested in the coarse scale part of the solution. The goal of our 
research here is to find an efficient way to compute it. 

Different frameworks for numerical multiscale methods have recently been 
proposed. The technique we will consider we will follow here is the heteroge- 
neous multiscale method (HMM) [51 H] . These multiscale methods couple 
simulations on macro- and microscales and compute the coarse scale solution 
efficiently In this paper we use HMM for wave propagation and we build on [9 
where we described a HMM multiscale method which captured the coarse scale 
behavior of ([!]) for finite time. The main aim here is to show that the same 
HMM methodology works also for long time, where new macroscale phenomena 
occurs. 

For finite time and periodic velocity coefficients converges to the solution of 
the homogenized equations, [5], 

U tt - V • (AVu) = 0, Qx[0,T], 
\u(x,0) = f(x), u t (x,0) =g(x), VzeQ, 

for fixed T independent of e, A E (x) = A(x,x/e) and where A(x,y) is periodic 
in y. The macroscale algorithm in HMM is inspired by ^ but it only does not 
use the explicit form and only uses the flux form, [9], 

( Utt -V.F = 0, nx[0,T% 
\u(x,0) = f(x), u t (x,0) = g(x), Vxefl. 

The challenge is to see if the same HMM approximation is robust such that 
it also is applicable for very long time, 0(e~ 2 ), where 0(1) dispersive effects 
appear. These effects are not given by |2]). The dispersive property is accurately 
captured if care is taken to have enough accuracy in each component of HMM. 

The numerical long time solution compares very well with analytic long time 
effective equation given by Santosa and Symes, [TS] . In the one-dimensional case, 
when A £ (x) = A{x/e) and A periodic, their long time effective equation will be 
of the form 

utt - Au xx - j3e 2 u xxxx = 0, (4) 

where A is the same coefficient as in |2]) and (3 is a functional of A which can 
be difficult to determine analytically and numerically. 

The rest of the paper is organized as follows. Section [2] contains a descrip- 
tion of the numerical algorithm including the added accuracy compared to the 
technique given in [S]. Analysis of the method is presented in Section [3] and the 
numerical results are contained in the following section. The conclusion are in 
the final Section [5l 
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2 HMM for the Wave Equation. 



We will now briefly describe how HMM can be applied to Q. We assume that 
there exists two models, a micro model h(u 6 , d e ) = describing the full problem, 
where u e is the quantity of interest and d e is the problem data (i.e. initial 
conditions, boundary conditions, . . . ), and a coarse macro model H(u, d) = 0, 
with solution u and data d. The micro model is accurate but is expensive to 
compute by traditional methods; in our case this model is The macro 

model gives a coarse scale solution u, assumed to be a good approximation 
of the microscale solution if and is less expensive to compute. The model is 
however incomplete in some sense and requires additional data. In our case we 
use 

u tt - V • F = 

with the flux F unknown. This is inspired by the form of the homogenized 
equation A key idea in the HMM is to provide the missing data in the 
macro model using a local solution to the micro model. Here ([l]) is solved 
locally on a small domain with size proportional to e and F is given as an 
average of the resulting microscopic flux A e \7u E . The initial data and boundary 
conditions (d 6 ) for this computation is constrained by the macroscale solution 
u. 

It should be noted that even if our numerical methods use ideas from ho- 
mogenization theory they do not solve any effective (e.g. homogenization or 
effective medium) equation directly. The goal is to develop computational tech- 
niques that can be used when there is no fully known macroscopic PDE. 

We will now describe a HMM method for the wave equation ([I]) which will 
give useful solutions in two regimes. The first regime is when T is fixed and 
independent of e. The other regime is when T — > oo as e — > 0. We will call 
this the long time regime and the problem itself a long time wave propagation 
problem. In this paper we will consider the one dimensional wave propagation 
problem. Many of the results can be shown to hold in a d-dimensional setting. 
We have in previous works shown higher dimensional results, both theoretical 
and numerical. We demonstrated in the previous papers [TU] that our HMM 
captures the same solution as homogenization (when applicable). In this paper 
we will primarily investigate how our HMM method handles the long time prob- 
lem. The microscopic variations in the medium introduces dispersive effects in 
the macroscopic behavior of the solution which becomes notable after long time. 
Our goal is to show that our HMM method can capture the dispersion with less 
computational cost than just resolving the full equation. 

Let us illustrate the dispersive effects by a numerical example. We consider 
a one-dimensional example where we solved with Q, = [0, 1], e = 0.01 and, 

( A £ {x) = A(x/e), A(y) = l.l + S3n2iry, Viefi, 
ju^O^e-^+e- 100 ^ 2 , ut(x,0)=0, VieO, 

for a long time, T £ — <D(e~ 2 ). Since we have periodic boundary conditions, 
the solution to the corresponding homogenized equation will be periodic in time 
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with period l/VA = 1.47722. We will compute the solution to for 1, 10 
and 100 periods. We see in Fig. [I] that after 100 periods there is an 0(1) error 
between the true solution u E and the homogenized solution u which thus fails 
to capture the dispersive behavior of the solution after long time. 
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Figure 1: Finite difference computation of ^ at T = 1.47722, T = 14.7722 
and T — 147.722 (1, 10 and 100 periods of the homogenized solution) and the 
corresponding homogenized solution (circles). As we can see the homogenized 
solution does not capture the dispersive effects that occur. 



2.1 Description of the Method 

The HMM method we suggest here is described in three separate steps. We 
follow the same strategy as in [T] for parabolic equations and in [2] for the one- 
dimensional advection equation. See [8], [5] and [5] for additional details and 
proofs. In step one we give the macroscopic PDE (i.e. the form H(u, d) = 0) and 
a corresponding numerical discretization. In step two we describe the microscale 
problem (microproblem) . The initial data for the microproblem is based on local 
macroscopic data. Finally, in step three we describe how we approximate F from 
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the computed microproblcm by taking a weighted average of its solution. 
We focus in this paper on the ID version of ([TJ, 

(i4 t -d x {A e ul)=0, ttx[0,T e ], () 
\u £ (x,0) = f(x), ut(x,0) = g(x), feed. 

where C K is a interval of the form [a, b] and u £ (x,t) is f2-periodic in x. 
Remark 1. All integrals are taken overM. unless explicitly stated otherwise. 

Step 1: Macro model and discretization We suppose there exists a 
macroscale PDE of the form, 

1^-^ = 0, Yx[0,T% 

\u(x,0) = f(x), u t (x,0) = g(x), Vxey, 

where u is ^-periodic; F is a function of x, u and higher derivatives of u. 
We will use this assumption throughout the whole paper. The assumption on 
([7| is that u ~ u £ when e is small. In the method we suppose that F — 
F{x,u,u Xl u xx ,u xxx ). In the clean homogenization case we would have F = 
Au x , and in the effective equation Q F — Au x + /3e 2 u XXX7 but we will not 
assume knowledge of the exact form of a homogenized equation or any other 
effective equation. 

We discretize Q using central differences with time step K and spatial grid 
size H in all directions, 

K 2 / 

rjn+l OTT n 7V n ~l J- I — TP n _l_ _ 07W n _l_ P" 

u m — Lu rn u m 2AH \ m + 3 / 2 m + 1 / 2 r m-l/2 ' r m-3/2 

(8) 

where F n j is F evaluated at and so forth. The scheme is second 

m— f m ~ 2 

order in K and fourth order in H . We choose to use a fourth order scheme in 
space since it showed to have better dispersive properties to allow us to avoid 
some of the numerical dispersion that pollutes the numerical solution. 

We will show how F is computed on the grid in the description of the micro 
problem. The implementation details are described in j!2| . 

Step 2: Micro problem The evaluation of F in each grid point is done by 
solving a micro problem to fill in the missing data in the macro model. Given 
the parameters x' and an initial data Q(x), we solve 

(vt t -d x {A%x + x>)v x )=0, r e x[-r,r], 

\ v £ {x,0) = Q{x + x'), v%(x,0)=0, VieT, U 

where v £ — v £ (x, 0) is y £ -periodic and Q depends on the macroscopic state {w^}, 
assuming x' — x m _ l / 2 , it is typically a third order polynomial interpolating 
iC_ 2 , ■ ■ ■ , < + i in x m - 2 , • ■ • , x m+ i, where u™ w U (x m ,t n ) in Q. Note that 
we have made a change of variables to center Y e around x' via x — x' i— > x. We 
keep the sides of the micro box Y £ of order e. We will discuss the choice of Y £ 
and Q below. 
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Step 3: Reconstruction step After we have solved for v £ for all Y e x [— r, r] 

we approximate F by a weighted average of the microscopic flux A e (x + x')v x 
over [ — 77, 77] x [— r, r] where [—77,17] C Y £ . We choose 77 sufficiently small for 
information not to have propagated into the region [—77, 77] from the boundary 
of the micro box Y E . More precisely, we consider smooth averaging kernels K„ 



compactly supported in [77,77], see below in Section 2.3.1 We then approximate 
FnF HMM (x',Q) = JJ K T {t)K v (x)f(x 1 t)dxdt, (10) 
where f e {x,t) = A E (x + x')v x (x,t) and v e solves Q. 
2.2 Motivation of the Method 

In this section we will give motivation to the steps and procedures of the HMM 
method for the wave equation, in the previous section. To do this we dehne the 
coarse scale solution u(t, x) as a local average in time and space of the fine scale 
solution, 

u(x,t) = ()Cu £ )(x,t), (11) 
where K is a local averaging operator defined as, 

(Kv)(x,t) = (K v * (K T *v))(x,t) = [[ K v (OK T (s)v(x-£,t-s)dtds. (12) 



The functions K v and K T are smooth, compactly supported, averaging kernels 
described below in Section [2.3.1| The operator K. commutes with d t and d x on 
smooth functions, since 



d t JCv = d t II K„(£)K r (a)v(x-£,t-8)d£d8 

K v {£)K T (s)v t {x - £, t - s) d£ds = JCv t . (13) 



In the same way we have d x Kv = K,v x . 

The coarse scale solution u will then satisfy a PDE of the form u u — d x F = 0, 
which can be seen as follows. We apply K. to both sides of the wave equation 

Ku\ t - Kd x (A e u e x ) = 0. (14) 

By applying the commutation property on <9 t /C twice and once for d x JC, we arrive 
at 

d tt {Ku £ )-d x {K{A s u x )) = Q, (15) 
which is the PDE for u that we seek, 

(u tt -d x F = 0, F = K,{A E u x ), 

I u(x,0) = (ICu £ )(x,0), ut(x,0) = (Kut)(x,0). [ ' 
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Note that the flux F is precisely the same type of flux used in Step 3 above. 
Hence, if u £ is the exact fine scale solution, the corresponding exact coarse scale 
solution must satisfy this PDE. To make the HMM procedure consistent, the 
same PDE is used for the micro problem. Additionally, the initial condition 
should be consistent, 

u(x,0) = (!Cu e )(x,0). (17) 

In Step 2 of the HMM procedure we are given the desired macroscopic initial 
data, i.e. u(x, 0). The initial condition Q(x) for the microscopic solution u e 



must then be chosen such that ( 17) holds. To leading order in e one can simply 
take Q(x) = u(x, 0), but to have higher order accuracy Q(x) must be modified. 
On the other hand, it turns out that the initial data for uf does not matter. We 
will discuss this consistency problem below in Section |2.3.2| 

2.3 Elements of the Method 

Now we will describe two important factors in the HMM method: The choice of 
kernel and how to obtain consistent microscopic initial data for a given macro- 
scopic data. We start with a discussion about kernels. 

2.3.1 Kernels 

In this section we introduce smooth compactly supported averaging kernels used 
to accurately compute averages of periodic and almost periodic functions. We 
will give a short motivation why we choose to work with those instead of just 
taking the usual mean value of / = t^t J n f(x) dx. We will also talk about how 
some numerical difficulties can be avoided by keeping the kernel on a factorized 
form. 

Suppose v e (x, t) is periodic function in x and t, with mean value v, and we 
wish to compute v with a numerical method. If we know the period of v e in 
both x and t we can compute v with spectral accuracy by the trapezoidal rule. 
The situation becomes harder when we do not know the exact period in one of 
or both of x and t. Suppose we have the values of v E available over some domain 
D = [—77, 77] x [ — t, t] where both rj and r are larger than but of the same order 
as e. The normal mean value then gives a first order error in e/rj and e/t, 



— r r « £ (£, S )d£d S 

VT J-r J-r, 



< C ~ + ~ ■ (18) 



To improve this error we follow the steps in |11| and introduce kernel functions 
of the following form. First, let K £ C%(M) satisfy 

K(x)x j dx =\ 1, 3 ~ °' supp-K" C [-1, 1]. (19) 

(0, j = l,...,p, 

We let K p,q denote the space of all such functions. In the computations we will 
scale K and write 

K v {x) = -k( X 
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This function inherits the integral properties of K in ( 19 1 , but supp K n c [ — 77, 77] . 
The improvement from normal averaging that one can achieve with K n is shown 
in the next theorem, which is a simple adaptation of some results derived in 

Theorem 1. Let f(t, s) be a 1-periodic function in the second variable whose 
derivative d r t f(t, s) is bounded for r = 0, . . . ,p. Then for f £ (t) :— f(t, t/e) and 
K G K p ' q we have 



K v (T)f e (t + T)dT~ / f(t, S )ds 



< drf + C 2 



where C\ and C'2 depends on the derivatives of f and K but not on rj or e. 
Moreover, if f(t,s) = a(t)b(s) and b is 1-periodic, then 



K n (r)f e (t + r)d7 



f(t, s)ds 



< C max ||a (r) | 

0<r<q 



L°°[t-r],t+rj\ 



(20) 



Hence, the convergence rate in ( 18 1 with respect to r\ and r can be improved 
significantly, 



K v {OK T {s)v%i,s)^ds 



< C 



(21) 



In our method we have chosen to use functions in K p,? which are polynomials 
and symmetric, K(x) — K(—x). Since the first q derivatives of K must vanish 
at x = ±1, those polynomials can always be factorized as 

K{x) = {l-x 2 ) q+1 P(x), (22) 

where P is a degree p — 1 polynomial. Finding the coefficients for P(x) and 
subsequently evaluating K(x) using the form (22 1 is much more stable numer- 
ically, than directly finding the coefficients for K, in particular for large p and 
q. This is because the coefficients of K then become very large. A code which 
computes the P coefficients for any p, q can be found in |12j . An example why 
it is important for the numerical accuracy to factorize the polynomial is shown 
in Table □ 



2.3.2 Consistency of the Microproblem Initial Data 

An important aspect in our HMM method is that the initial data for the mi- 
cro problem is consistent with the current macroscopic state. In practice we 
approximate the macroscopic state with a third order polynomial U(x) which 
interpolates the macroscopic grid values • • • > u m+i m x m-2, ■ ■ ■ , %m+i- 

We should then choose a polynomial Q{x) as initial data to ^ such that 

U{x) = (K,v E )(0,x). 

To simplify the discussion we introduce the operator M. which maps the initial 
data of (JoJ) to (K.v E )(0, x). We hence seek Q such that 

U{x) = (MQ)(x). 
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N ^^^^ 


14 terms 
, s 

7475.95 • ■ • x 28 + ■ ■ ■ 


(1 - x 2 Y +l P{x) 


10 


6.0572e-04 


6.0572e-04 


20 


1.5465e-07 


1.5467e-07 


40 


2.2333e-10 


2.6219e-10 


80 


3.6733e-ll 


8.3822e-14 


160 


4.8612c-ll 


1.1102e-16 



Table 1: Here we have plotted the relative error as a function of N (grid inter- 
vals) for two ways of evaluating the integral of the same polynomial K £ K 9 ' 9 . 
The polynomial evalutions themselves are done with Horner's method (MAT- 
LAB builtin function polyval). We clearly see that there is a great benefit in 
using the factorized K from a numerical perspective. The 28 degree polynomial 
K has only even powers of x l , but has very large coefficients. The large coeffi- 
cents and the number of terms makes it difficult to evaulate K accuractly. The 
expected convergence rate is spectral. 



We will not be able to find a Q where this is satified exactly. However, in this 
section we will show how to find Q given U such this equality is satisfied to high 
order in s. 

When the coefficient A £ is precisely periodic and Q is a third order polyno- 
mial, we show below in Section [3. 3| that, under some simplifying assumptions, 

(MQ){x) = Q(x) + e 2 1 Q"{x) + G(e 3 ), (23) 

where 7 is a constant. Hence, if we take Q = U we make a 0{e 2 ) error. For the 
long time problem we will need better accuracy. 

We assume that the initial data is of the form Q(x) — cq + c\x + C2X 2 + C32 
and make that ansatz that 

(MQ)(x) ~ U(x) := c + cix + c 2 x 2 + c 3 x 3 . 



[cq ci c 2 c 3 ] , c = [c ci c 2 63] and x T = [l 



x x 



Then Q(x) = x T c and U(x) = x T c. We define as the unit vector with a one 
in the ith component (counting from zero) and / as the 4x4 identity matrix. 
For consistency we want to figure out which c that gives us a particular c and 
then compute the flux F for that c. As an example, assume that Q{x) = x 3 , 
which would mean that c = e 3 . We cannot use this as initial data to obtain 



the correct flux as it will be inconsistent (to order e 2 ), as seen from (23) and as 
described briefly in 10J. Instead we wish to figure out which initial data x T c 
that gives us a x T c such that c = e 3 and use that x T c as initial data. We will 
show that it is not necessary to recompute any micro problems once we have 
solved just four independent micro problems. 

Since Ai is a linear operator we can represent it by a matrix M acting 
between the coefficients c and c, 

c = Mc. (24) 
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We now note that we can find the elements of M from solving four separate 
microproblems as follows. For i = 0, . . . , 3 we take Qi(x) Solving 
the corresponding microproblcm gives Ui := AAQi, which we can evaluate in 
four points. We fit a third order polynomial to these points and take this as an 
approximation of Ui(x) —: x T Ci. Then, is the i-th column of M, 



Me, 



m ,i 

m 2 .i 

m 3 ,i 



Once we have found M we can solve (24) to obtain c from c. We call M 
the correction matrix. Based on the scaling in ( 23 ) we can see that M is a 



e -perturbation of the identity, M = I + e Mq for some Mq. 

Once we have computed the M matrix we can also compute the correct flux, 
without solving more micro problems. For each of the four problems solved 
we record the corresponding flux and denote it by /j. Moreover, we set f T = 
[ fo fi fi J3] ■ Since the flux computation is also linear in the input c the 
computed flux from Q(x) is simply 



HQ) 



Thus, given U = x c the corrected flux is 



F(U) = f c, f-^M-'f. 



(25) 



Let N be the number of grid points on the shifted grid x m _i (1 < m < N). 

We determine the coefficients in the 4x4 matrices correction matrices 
with the following algorithm 
for to = 1 — >• N do 

for i = 1 — > 4 do 



Initial data Q(x) <— x % , i.e. 



1 and other c, = 



Solve micro problem ([9| v £ for all (x,t) € [ — to rj\ x [— r, r] 
Compute Yk K,(u £ )(Ak, 0) for some distinct around x = 
Find a third degree polynomial that fits (Afc, Yfe), which gives c. 



M (m) : 

end for 



Cj-l, 



1 < i < 4 



end for 

So far we have only been concerned with the initial data for v e . In ^ we 
also need initial data for uf . However, as the following derivation shows, the v\ 
data does not affect the flux when we use symmetric kernels. We consider the 
wave propagation problem of the form 



vl t -d x (A s v s x )^0, 
v s (x,0) = f(x), v s t (x,0) 



(26) 
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We can split the solution into two parts V s = vf + tj| where vf solves ( 26 ) with 
the initial conditions 



vl(x) = f(x), d t vl(x,0)=0, 



and tj| solves ( 26 ) with the initial conditions 

«!(*)= 0, d t v s 2 (x,0) = g(x). 

The solution u| will be anti-symmetric in time, i.e. v^—t^x) = — 7j|(f, x), since 
A £ is time- independent. The same holds for A £ d x V2- Then, because our kernels 
are symmetric, K,(A e d x V2.)(0, x) = and 

F(0, 0) = K(A E d x v e )(0, 0) = IC(A £ d x vl)(0, 0), 

independent ol .9(2;). In the computations we therefore simply take g = 0. 

2.4 Changes of the HMM Method for Finite Time for Ac- 
curate Long Time Wave Propagation 

We have have made three important changes to the HMM method for finite 
time [5] to capture the effective solution for long time wave propagation: 

1. The form of the macroscopic PDE Q is the same, u u — d x F — 0, but the 
flux F = Au x + (3e 2 u xxx contains a third derivative of the macroscopic 
solution. To capture the additional information we need to pass more 
information to the micro problem. The initial data is chosen as a third 
order polynomial, instead of a linear polynomial as for the HMM method 
for finite time. 

2. We need to make sure that the initial data in the micro problem is con- 
sistent to a very high accuracy with the macroscopic solution. As we 
described in Section |2.3.2[ the initial data to the micro problem needs to 
be consistent with the macroscopic data to accurately capture the disper- 
sive effects which are important for long time wave propagation. 

3. Since we need to accurately represent also the second term in the flux, 
j3e 2 u xxx , the error in the flux computation (10) must be smaller than 



C(e 2 ). By Theorem [T] the error £ in the flux computation is proportional 
to, 

£~ V p + + {e/TY + (e/r,y. 

Thus, we need to choose p, q, r\ and 77 such that £ <C e 2 , asymptotically. 
This means that the micro box must be larger than e. For instance, if 
p = q > 2 we can take take a < 1 — 2/q. Then 

^ 1 T~s a => £ = O (e p - 2 < 0(e 2 ). 

Recall that in the finite time case we always take 77, r ~ e, This hence 
explains why we need to have more accurate kernels or bigger micro boxes 
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in the long time case. In order to maintain a low computational cost 
we usually want small micro box, which can be obtained by using a very 
regular kernel, i.e. large q. 

Remark 2. We use the same numerical scheme ^ to solve the micro problem 
([9|, where is computed as, 



F v 



A*(x % 



24h 



h3/2 ' 



27vl 



.+ 1/2 



-1/2 



1-3/2 



(27) 



3 Theory 

In this section we will analyze the flux used in HMM. In particular we want to 
compare it with the corresponding flux for the effective equation. We will focus 
on the case where we have a good understanding of this equation, namely when 
the coefficient A e [x) is precisely periodic. We thus consider the micro problem 
centered at x = extended to the whole space, 

{ " V ' ' . (28) 

| v e (0,x) = Q(x), u|(0,a;) = 0, i6l , 

where A(y) is 1-periodic. In [9] we studied the finite time case where Q(x) was 
a first order polynomial. We could then prove that 

F H MM(Q,a)=F HOM (Q) + 0{a q ), -Fhom(Q) := AS/Q, (29) 

where we have defined 

e 

a = -. 
V 

This means that when the microbox size r\ becomes large in comparison to e, 
hence when a — > 0, the HMM procedure generates results close to a direct 
discretization of the homogenized equation ^ . 

The restriction to linear initial data made the analysis significantly easier, 
since then v £ (t, x)—v e (0, x) is periodic for all t and the solution could be analyzed 
via a rather straightforward spectral decomposition. For the long time case we 
need to use initial data Q(x) which is at least a cubic polynomial. To analyze 
this case we need to use the full Bloch theory for wave propagation in a periodic 
medium. This was also used by Santosa and Symes in [T5] to describe the 
effective properties for periodic medium over long time. We start by reviewing 
this theory. 



3.1 Long Time Homogenization 

The theory for long time homogenization formally extends the validity of the 
effective model ^ from T = 0(1) up to time T = 0(e~ 2 ). It does this by 
adding more terms to the effective equation. The basic model was presented in 
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We start with some notation and definitions for Bloch waves. We let 
and ipm be the eigenvalues and eigenfunctions of the shifted cell (eigenvalue) 
problem [3l pp. 614], 

f - (V„ + ik) A(y) (V y + ik) $(y, k) = u 2 {k)^{ Vl k), Y x Y, 
1 ip(y, k) is F-periodic in y, 

where Y — [0, l] d and k G R d . Let 4>m(x, k) be the e-scaled Bloch- waves, 
<f> m (x, k) = ip m (x/e, ek) exp(ifcz;), 

which satisfies 



We will now consider the general solution to ( 28 1 when we assume that the 
initial data Q(x) is band limited with smallest wavelength larger than e, 

Q(x) = J Q(k)e ixk dk, supp Q C [-tt/e.tt/e]. (30) 

We let V m and q m be the parts of v e and Q that project on <p m , 

V m (t,k)= v e (t,x)<f)* rn (x,k)dx, q m (k) = Q(x)4>* m (x, k) dx. 



Then 

OO „ 

" " id*. 



v £ (t,x) = V / V m (t,k)<t)m(x,k)i 



Upon inserting this in ( 28 1 we see that V m will satisfy the ODE 

dttVm + \u? m {ek)V m = 0, V m (Q, k) = q m {k), d t V m (0, k) = 0. 



By solving this we can then write down the exact solution to ( 28 ) as 

OO r, OO 

v e (t, x)=^2 / q m (k)4> m (x,k)cos(uj m (ek)t/s)dk -. ^2 v m {t,x), (31) 
where 



v m (t,x)= / q m (k)4> m (x,k)cos(u} m (ek)t/e)dk. 
Jy/s 

The following theorem from [15] states that vq is the leading order contribution 
to v £ . In fact the sum of the terms with to > 1 is bounded by 0{e) in L 2 , 
uniformly in time. 
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Theorem 2 (Santosa & Symes). Suppose v € solves (28) with d = 3 and initial 
Q satisfying \3(ty. Then 



m — 1 



dx < Ce 2 



Here C is independent of e and t but depends on the H 2 -norm of the initial data 
f and the L°° -norm of A and dA. 

See [T5] for proof. The long time homogenized equation is an equation, 
approximately, satisfied by the leading term vq derived in |15j . In one dimension 
it is of the form 

u tt = Au xx + e 2 /3u xxxx . (32) 
The homogenized coefficients A and (3 are given as derivatives of 

n(fe) := uio(k) 2 , 

not to be confused with the spatial domain f2, evaluated at k = 0, 



A=io"(Q), 







~n""(o). 



(33) 



This model is in principle valid up to time t = 0(s~ 2 ). However, the equation 
is not well-posed since the sign of /3 is wrong. When the grid is refined the 
numerical solution becomes unstable. Still, solving this numerically on a coarse 
grid gives solutions that agrees well with the long time behaviour of (28) as 
demonstrated in |T5j [10] ■ See [10] for further discussion and analysis of this 
issue. 

Remark 3. There is another formulation of the effective equation which is 
also valid up to T = 0{e~ 2 ). This equation is well-posed but the higher order 
correction term involves time derivatives. It is of the form 



i 20 
u tt = Au xx + e -j 



Hxxtt • 



It was recently analyzed in 
also for the long time case. 



where convergence was established rigorously 



3.2 The Flux 



As seen in ( 32 ) and ( 33 1 the homogenized flux is 

fhWQ) = AQ' + (3e 2 Q"' = ^"(O)Q' - £ 2 ift""(0)Q" 



(34) 



We will here analyze the expression for the corresponding HMM flux in the one- 
dimensional case. This flux is computed using the averaging kernels from Section 



2.3.1 Let v e be a solution to the HMM micro problem (28|. Henceforth we 
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drop the explicit dependence on initial data Q in the notation for -Fhmm(Q, °)- 
Then, 



m=0 

1-1/26 



^HMM(a) := / / K v (x)K 11 (t)a(x/e)d x v £ (t, x)dxdt 

oo r r 

= J2 I K v (x)K r ,{t)a(x/s)x 

q m (k)d x ip m (pefc) e lkx cos ^ m ^ k ^ dk dxdt. (35) 

The goal is to study how -FhmmO*) behaves as a — ¥ 0, hence, when the microbox 
size rj becomes large in comparison to s. (We always assume a < 1.) We first 
note that for any r(x), the expression 

/KJx^rix I e) dx = [ -K(x/i])r(x/e) dx = — [ K(xe/rj)r(x) dx 
J V V J 

= ajK(xa)r(x)dx, 
is a function of the ratio a = e/rj, not of r\ or e individually. We therefore set 

b m (k) = [ ^* m (x,k)dx, s m (fc,a) = / K v (t) cos(uj m (k)t / e) dt. 

J -1/2 J 

Moreover, let 

w m (k,a)=e J K, 1 (x)a(x/e)d x ip m (x/e,k)e tkx/e dx, 
and finally define 

oo 

W(k,a) = b m (k)s m (k, a)w m (k,a). (36) 

m=0 

We can then prove the following theorem which, under a smoothness assumption 
on W, shows that when initial data is a polynomial the HMM flux (scaled by e) 
is always a polynomial in ed x of the same order, whose coefficients only depend 
on a, not directly on e. 

Theorem 3. Suppose a G C 1 , W(-,a) € L 2 {— ir,ir) and |fffl| ) holds. Then 
FnMu(a) = - [ Q(k)W(ek,a)dk. 



e _ 

Moreover, if Q is a polynomial of degree r and W(-,a) G C r {— tt,tt), 

Fhmm = £^ 1 a o (a)Q(0) + oi(a)Q'(0) + ea 2 (a)Q"(0) + • • ■ + s r - l a r (a)d r x Q{Q). 

where 



a r (a)=' > —dlW(0,a). (37) 
r! 
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Proof. We prove this theorem in several steps, beginning with a lemma. 
Lemma 1. Suppose the support of Q lies in [—ir/e,ir/s\. Then 

q m (k) = Q(k)b m (sk), 

and in particular supp q m C supp Q. 

Proof. Since ip m (x,k) is 1-periodic in x and belongs to L 2 for each fixed k we 
can expand it in a Fourier series, 



i; m (x, fc) = £ b* m Jk)e-* 2 ™ x , b m , n {k) = f' 2 r m (x, k)e~^ nx dx. 

„ J-l/2 



,1/2 
1/2 

Then 

-ikx 



q m (k) = J Q(x)ip* m (x/e,ek)e lKX dx 

= EW(^) / Q(x)e-« k - 2 ™^ x dx = J2 b m,n(ck)Q(k -2i.nl 'e)). 



But the support of Q is in [— n/e, n/e] and when n^O, 



27m 
k 

e 



. , , 27T 7T 2-7T 7T 

> - k + —> — + — = -, 

e see 



so the only contribution to the sum is when n = 0. Hence, 

q m (k) = b mfi (ek)Q(k), 
and the result follows. □ 
From the result of this lemma we then get 

J2 K ri {x)K v (t)a(x/e)q m {k)x 
d x ipm{x/£ 7 ek) exp(ikx) cos{uj m (ek)t/e) dx dtdk 
Q(k) b m (ek)s m (ek,a)x 



tt/s 



m=0 



J K 11 (x)a(x/e)d x ip m {x/e 1 ek) exp(iekx/e) dx dk 
: ~ / Q(k) V" b m (ek)s m (ek,a)w m {ek,a)dk. 



This gives the result in the first part of the theorem. To show the second 
part of the theorem, we note that since the mapping Q \-t Fhmm is linear, we 
only need to prove the result for Q(x) = x r . Polynomials are not in 1? so we 
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first regularize Q. Let M £ £ C£°(R) have the properties that M e (k) = for 
|fc| > n/e, M s (0) = 1 and M £ W (0) = for £ = 1, . . . , r. Then define 



Q„(aO= // Afe^x/^^yV^-^dydfc, 



-a, a 



By construction this function is in L 2 and band limited to [— 7r/e, tt/e] for all val- 
ues of a. We let Q a denote its Fourier transform and F a the flux corresponding 
to initial data Q a . Then 



Q a (k) = f Me {k) — D a {k), D a (k) 



-iky 



dy. 



and 



F a = 



] a (k)W(ek,a) dk 



D^(k)M e (k)W(ek,a) dk 



D a (k)—M e (k)W ( - r \ek,a)dk. 
dk r 



This is valid for all a. We can then take a — > oo. Since D a (k) — ^ S(k) and the 
derivatives of M F vanishes at k = 0, 



15m F a = 

a->oo £ clfc r 



(fe)' 



T^ (r) (0,a). 



k=0 



Moreover, 

lim Q Q (a;)= Um / Q a (k)e ikx dk = i r lim / M e (k)Di r) (k)e tkx dk 



= (-i)''(ia;) r =x r . 
This proves the theorem. 



fc=0 



□ 



3.2.1 Properties of VT(fc, a) 

In this section we prove a theorem about W(fc, a). What we would like to show 
is that the coefficients ( 37 1 in Theorem [3] satisfy do (a) = and 

a r (a) = ~W0,a) = £L ^n(k)B(k) 



where we have defined 



(r + iy. dk r+1 



k=0 



-0(e- r+1 a q ), r > 1, 
(38) 



:= \b (k)\ 2 



In other words, we would have 
i d r+1 



d r k W(0,a) 



r + 1 dk r + 



Y n(k)B(k) 



+ 0(e^ +1 a 9 ), 



r > 1. 



fc=0 
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This would be consistent with our numerical experiments and with ( 34 ) for r < 2 
when a — > 0. We can see this in the following way. Recalling that wnfO) = wc 
note that fi(0) and f2( p )(0) = for p odd. Moreover, from Lemma |2| below we 
also have B'(0) = and B(0) = 1. Hence, 



ai ( a ) = -(n"(0)B(Q) + 2fi / (0)B / (0)+n(0)B // (0))+ 



jfi"(0)+O(a«), 



and 



ea 2 (a) 



3! 



(fi'"B + 3fi"B' + 3n'fl" + OB'") + C(a 9 ) = 0(a 9 ). 



For r = 3 have already noted that initial data must be chosen more carefully to 
have consistency. The value of 03(a) must therefore be different from what is 
given in (34). Indeed, 



e 2 a 3 (a) 



4! 



4! 



(Q""B + 4n"'B' + 6n"B" + iQ'B'" + VLB"") + 0{a q ) 
(fi""(0) + 6J2"(0)B"(0)) + 0(a q ). 



Hence, in total, using Theorem [3] and the expressions for a r , when Q is a third 
order polynomial, 

FhmmW = ^O"(0)Q'(0) - ^ [Q""(0) + 6O"(0)B"(0)] Q"'(0) + O(a'). 



The difference with ( 34 ) is in the second term, which comes in because the initial 
data is not consistent on this level with the macro data. The difference from 
the case of linear initial data, which was analyzed in [5] , is clearly seen; when 
Q'" = the flux is 0(a q ) away from the homogenized flux as in (29). However, 
we will see below in Section 13.31 that if the initial data is made consistent then 
the result agrees precisely with (34 1 if 03 (a) is given by (38 1 . 



Unfortunately, we are, not able to prove the full result in ( 38 1 at this point 
and it stands as a conjecture. Nevertheless, we have some partial results sum- 
marized in the theorem below. In particular we can show that (38) is true for 
r = 0, 1. For r > 2 it is true for the first term in the infinite series that defines 
W(k,a) in @, 

Wa(k,a) := bQ(k)so(k, a)wo(k, a). 



This term corresponds to the leading order part of the solution vq in (31). Wc 
have 

Theorem 4. The function W(-,a) belongs to L 2 (—tt,7t) D L°°(— it, it) with 
norms bounded independently of a < a®. When K G W ,q and W(-,a) £ 
C r {— 7T, 71") we have for r = 0, 1. 



W(0,a) = 0, VF (1) (0,a) 



i_ _d 
2 dfc 



2 n(k)B(k) 



+ 0(a q ). 



fe=0 
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and for r > 2, 

+ 0(e- r+1 a q ). (39) 



Wo (r) (o.«) = ^^n(*)s(*) 



fe=0 



To prove this we first need to establish some lemmas about the coefficients 
b m , s m and w m , starting with the b m coefficients. 

Lemma 2. We have 

oo oo 

^|6 m (fc)| 2 -l, b m (0)=5 m , J2 b 'm(°)'f , m(x,0) = -d k Mx,0)- 

m—0 m— 

Moreover, for the squared quantity B(k) = |o (fc)| 2 , 

B(0) = 1, B'(0) = 0. 

Proof. By definition, 

oo 

1 = b m(k)ip m (x,k), 

m=0 

since {tp m } is an orthonormal L 2 -basis for each fixed k. Taking the L 2 norm over 
[—1/2, 1/2] we obtain the first result. The second result is true since ipo{x,0) = 
1, and then 



,1/2 ,1/2 

b m (Q) = tp^ l (x,0)dx = ipl l (x,0)ip {x,0)dx = d m . 

J-l/2 J-l/2 

Together with, 



= d k ^2 b m (k)tp m (x,k) = ^2 b' m (k)i> m (x,k) + ^ b m (k)d k ip m (x,k), 

m= m—0 m—0 

the third result follows. That B(0) = 1 is true since 6 m (0) = <5 TO . For the last 
statement we note that for all k, 

= J^l = y ipo(x, k)ipo(x, k)dx = 23ft J tp (x,k)d k ipo(x,k)dx. 

In particular, for k = we get = 23?6 (0) since i/'o(a ; 1 0) = 1. We finally note 
that B'(0) = 23?6 (0)6 (0). □ 

Lemma 3. Suppose K e W' q . For s m (k,a) we have when m = 



k=0 



1, r = 0, 

0, l<r<p, 



and /or m > 1, 

|s m (0,a)| < Ca 9 , m>l, 
w/iere £/ie constant C is independent of a and m. 
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Proof. Since w (0) = 0, 



s o (0,a) = / K v (t)dt = l. 



For r > 1, by the moment conditions (19), 
d 



dk r 



so(k,a) 



fc=0 



/ K v {t) 



/ K v {t) 



it {it) 2 

Cl.r h C2,r o h • ' ' + c r r 

it (it) 2 

Cl,r h C 2 , r 5 !-•■•+ C r r 



(«)' 



exp(iw o (0)i/e) df 
dt = 0. 



Moreover, for m > 1, by Theorem [T] 



|s m (0,a)| 



^(t) cos(w m (0)t/e) dt 



< C 



27T£ V 

'm(o)??y 



<c" 



wi(0) 



□ 



Lemma 4. Suppose a(x) g C 1 arwi if € K p,? wrat/i g > 1. TTien 

sup \w m (k,a)\ < C, 

\k\<n 

with C independent of m and k. Let 

h(x, a) — — a 2 K 1 (ax)a(x) — aK(ax)a! (x). 
Then h € C c with supp h C [— 1/a, 1/a] and we can write 

ikx 



w m (k,a)= / h(x,a)ip m (x,k)e l x dx. 



Moreover, 



j \w m {k,a)\ 2 Ak= \h(x)\ 2 dx <C(l + a) 2 , 

where the constant C is independent of k and a. In addition, wo(0,a) = and 
for 1 < r < p, 



Proof. Let 



e 1 

:= — ed x K„(x)a(x/s) — — ~K (x I rj)a(x I e) K(x/rj)a (x/e). 

r\ l rj 
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Moreover, via intergration by parts, 

w m (k,a) = e j K v (x)a(x/e)d x ip m (x/e, k)e tkx/e dx (40) 

'h(x)^(x/e,k)e^ax. 
We have supp h C [— r], rj\ and \h\ < C(a + l)/?y. Therefore by Cauchy-Schwarz 



n 



\w m (k,a)\ 2 < / \h{x)\ 2 dxx / \ip m (x/s,k)\ 2 dx 

J Tj J Tj 

<C \a + l\ 2 /r] 2 dxxe \ip m (x,k)\ 2 dx 

Jr, Jl/a 

rj a 



This shows the boundedness. Next, from (40) it is clear that 

w m (k,a) = e I h{ex)ip m {x, k)e lkx dx = / h(x, a)ip m (x, k)e lkx dx. 



Hence, w m (k,a) are the Bloch coefficients for h(x,a). Therefore, by Parseval, 
J \w m (k,a)\ 2 dk = J \h(x,a)\ 2 dx = -J \h(x/a, a)\ 2 dx < C(l + a) 2 . 



m=0 



To prove the last result we let tp = tpo to simplify the notation. Denoting the 
binomial cofficients by c r i 1 we have 



dk 



: w Q (k,a) 



dk r 



K n (x)a(x/e)d x ip(x/£, k)e lkx/e dx 



r-l 



iJ2 c n / K v (x)a{x/e)d x ^ l \x/e 1 k)e lkx l 
en f K n (x)a(x/e)^ (x/e, k)e lkx ' 

1=0 

+ J2c rl f K v (x)a(x/e)ik^ e) (x/e,k)e tkx/E ( — 
1=0 J ^ £ 

r-l 

+ J2c re K ri {x)a{x/e)^ e) (x/e,k)e ikx/£ i{r - £) 



dx 



dx 



r-l 



dx 



1=0 



dx. 



Since ij)(x, 0) = 1, 

wi r) (0,a)= Y. c ^ / K v (x)a(x/e)^\x/e,0) 



dx 



1=0 



+ Y,cn / K„(x)a(x/e)il>W(x/e,0)i(r-£) ( — 



r-l-l 



dx. 
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We can now use ( 20 ) in Theorem [T] to get 



1/2 



w { o'(0,a) = / a{x){^ ) {x^) + ic rtr - 1 iP {r '' 1) {x,Q))Ax+ 0{e- r+1 a q ) 



1/2 

= Km — / afcXVfcfo fc) + k)) dx + 0{e- r+1 a q ) 

fe— >0 dfc r J— 1/2 

d r 1 f 1 / 2 

lim— — / (fix +i£>(x)(a E + ifc)V>(a;, fc) dx + 0(e _r+1 a 9 ) 



fe-vo dfc r ik 



1/2 



lim "°. (fc)2 xP(x,k)dx+ 0(e- r+1 a q ) 



fc-s-o dfc r ik 



-1/2 



fc->o dfc r ifc 
This concludes the proof. 



hm 4- ^° (fc)2 65(fc) + 0(e- r+1 ^). 



□ 



3.2.2 Proof of Theorem g] 

To show that W £ L 2 ( —ir, it) we note that by Lemma|2j Lemma[3]and Lemma[4] 



\W(k,a)\ 2 dk = 



2J b m (k)s m (k,a)w m (k,a) 

m—0 

/7T OO cxj 

sup|s m (fc,a)| 2 ^ |6 m (fc)| 2 ^ |w m (A:,a)| 2 dfc 
" m—0 m—0 

<C^2 \w m (k,a)\ 2 dk <C. 



m=0 ' 



We should also check that W £ L°°. Let 

OO 

g(x,k,a) = ^ b m (k)s m (k,a)ip m (x,k). 
Then with N = \l/a] 



m=0 



N 



W(k,a)= / h(x, a)g(x, k, a)e dx 

J-N 

and since \h\ < Ca and by the periodicity of g, 

/N nl/2 
\g(x,k,a)\dx<C / \g(x,k,a)\dx<C\\g(-,k,a)\\ L 2. 
-N J -1/2 
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We note that since {ip m } is an ortho-normal L 2 -basis for each fixed k and 
\s m (k,a)\ < C, 



/1/2 oo oo 

\g\ 2 dx = V \b m (k)\ 2 \s m (k,a)\ 2 < C 2 V \b m (k)\ 2 = C 2 . 
-1/2 ±Tn ±Tn 



Hence, |W^| < C uniformly in k. This shows that W € with a norm bounded 
uniformly in k and a. 

We have by the lemmas above, and since wo(0, a) = 0, 

oo 

W(0,a) = ^& m (0)s m (0 , a)-u; m (0, a) = 6 (0)s (0, a)w o (0, a) = 0. 

?TI = 

For 1 < r < p, with binomial coefficients c r g we get 



IT 1 

W / o (r) (0,a)= —b (k)s {k ia )wo(k,a) 



k=0 



= ^2c re s ( r e) (0) —j b Q {k)w (k,a) 

1=0 

d 



fe=0 



dk 



: b (k)w (k,a) 



k=0 



e=o 



E 



(r-l) 



(0) 



1=0 

■ d r fl(k)., ..... 



/A: 



fe=0 



+ 0(e- r+1 a 9 ) 



fc=0 



Since if 5 (/c) = fc/(fc) we have g {r) {k) = fc/M(fc) + rf^" 1 ^), 



d r n(k)B(k) 



dk r 



k 



1 



qr+l 



fc=0 



r + 1 dfc r+1 



n(k)B(k) 



k=0 



small. We set 



which proves (39 1. We should finally show that the difference — W^p is 



z(x, k, a) = [g(x, k, a) — bo(ot)so(k, a)ipo(x, k)] e 



ihx 



Then 



A' 



W(k, a) — Wo(k, a) = / h(x, a)z(x, k, a) dx 

J-N 

and since b m (0) — 5 rn by Lemma [2j 

oo 

z k (x,0,a) = ^2 b' m (0)s m (Q,a)ip m (x ) 0). 



in— 1 
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As before, 



JY 



1/2 



\W (1) {k,a) - W (1) (fc,a)| < Ca / \z k (x, k, a)\ dx < C / \z k (x,k,a)\dx 

J-N J -1/2 

<C\\z k (;0,a)\\ L 2. 

Furthermore, 



\ z kC 



(•MWh = E IC(0)| 2 MM)| a < Ca 2 i £ |C(0)f 



m— 1 



Hence, the size of the remaining sum is of the same order as the error 0(e r+1 a q ) 
in the leading term for r = 1. 

3.3 Consistency 



In this section we make a formal derivation to motivate the relationship ( 23 ) 



We therefore consider the exact solution for initial data Q(x) which was given 



above in (31 1. We assume that only the zeroth term vq is relevant and that the 
initial data Q is band limited as in Q. Then, by Lemma [I] and 



l/2e 



l/2s 



Q(k)b (sk)ipo(x/e, ek) exp(ifcx) cos(ujo(ek)t/e)dk. 



V £ (t,x) W Wo(*)3f) : 

This implies that 
where 

U(x,y) = / Q(k)bo(ek)sQ(ek, a)tpo(y, ek) exp(ikx)dk. 

Since A; is in a bounded set we can Taylor expand the terms in efc. From 
Lemma [3| 

s (£fc,a) = 1 + O ({ek) p+1 ) . 

Hence, 



P d l f 



Q(k)exp(ikx)dk+ 0(e p+1 ) 



E 

£=0 



(-is)* 



bo(k)ip (y, k) 



Q&{x)+ 0(e p+L ) 



fe=0 
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We can use Theorem [T] to deduce that 

,1/2 

K^x' -x)d r Mx' /e,0)Q {i) (x')&x' = / d r k ip (x^)dxQ {l) (x) + 0(a q ) 



-1/2 

= d r k b (0)*Q^(x) + 0(a q ). 

Therefore, 

p (— ' V 

(MQ)(x) w ^B ffl (0)Q (fl (i) + 0(e*a«) + 
bmce B(0) = 1 and B'(0) = by Lemma[2] we have, in particular when p > 2 
(Xg)(x) w Q(x) - yB"(0)Q"(!t) + 0(e 3 + a q ). 



This is the expression (231 with 7 = —B"(0)/2. In order to have an accuracy 
that is better than an 0(e 2 ) we need to correct Q(x) and to use a kernel such 
that a q = 0(e 3 ). More precisely, we should take 

Q(x) = Q(x) + ^-B"(0)Q"(x), 
and use this as initial data instead of Q. Then, if Q is a third order polynomial, 

(/Cu)(0,a;) « Q(i) - £ -B"{Q)Q"(x) + 0(e 3 ) 

Q(x) + yB"(0)Q"(x)j - ^2?"(0)Q"(:e) + 0(e 3 ) 
= Q(a;) + 0(£ 3 ). 

Let us now see what the implications are for the computation of the flux Fhmm- 
If Q(x) is a third order polynomial then so is Q(x) and we can use the result in 
the previous section. With initial data Q{x) we obtain 

1 O"" -I- fiO" R" 

Fhmm = Vq'(O) - £ 2 " + b» g q,„ (q) + 0(q9) 

1 / \ O"" 4- fiO" R" 

= in" r Q'(0) + y W(0)J - e 2 + 4 " Q w (0) + 0(« 9 ) 

= 2n"Q'(o)- £ 2 ^Q"'(o) + o( £ 3 ). 



This is thus consistent with (34 1 upto order 0(e) if ofl — 0(e ). 



4 Numerical Examples 

In this section we consider three long time wave propagation problems and a 
detailed convergence study of the flux correction. The three problems are of the 
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form, 



d x (A E u x )=0, [0,1] x [0,T% 



u £ (x,0) 



,-100(1-2: 



)2 , u e t (x,0) = 0, Vx€[0,l], 



(41) 



with 1-pcriodic boundary conditions in x and T E = 0(e~ 2 ). Note that the initial 
data has no e dependency. We will try three different A £ = A% parameters: 



1.1 



A\{x) 
A|(x) = 1.1 



2ttx 



cos 2nx + sin 



27rx\ 



A%{x) = (l.l + sin2vr^ (1.5 + 0.5(cos(2ttx - 1))) . 

The problems can be expected to have dispersive effects which are not described 
by the homogenized wave equation ([2| restated here for convenience, 



u tt - d x (Au x ) = 0, 



[0,1] x [0,T], 



u{x, 0) 



ox* +e -ioo(i-x) 2 j Ut (x,0)=0, Vx€[0,l], 



(42) 



where u is 1-periodic. We suppose there is an effective equation valid for T e 
0(e~ 2 ) of the form 



u tt - d x (A{x)u x + 0(x) 



= 0, 



[0,1] x [0,T% 



u{x, 0) 



^ +e -ioo(i-x ) 2 j ^05 = 0, Vie [0,1], 



(43) 



where A is the same A as in (42 1, but unfortunately as we described in relation to 
Q, (3 can be difficult to determine from A, both symbolically and numerically. 
We compute A(x) and j3(x) from a A £ of the form A e (x) — A(x, x/e) by freezing 
x and computing A(x) and j3(x) as if A e had only fast oscillations. We have 
used Maple to compute the f3 coefficients from A, [T2"] . 

We will consider four methods for each wave propagation problem: An exact 
(DNS) solution of (41 ) where we discretized the full problem with a finite differ- 
ence method; a discretization of (42 1, the corresponding homogenized equation 
(HOM); a discretization of (43), the corresponding effective equation (EFF); 
and finally a HMM solution. Throughout our examples we will use the follow- 
ing parameters: e — 0.03, rj = r = 20e and a polynomial kernel which is 19 
times continuously differentiable and has 19 zero moments. We will also use the 
notation, 

At £ e 1 

A:=— , p e := — , p:= 



Ax' 



Ax' 



Ax 



4.1 Wave Propagation Problem One 

We consider Af (x) = 1.1 + sin27ra;/e, which has only a fast scale. For the A\ at 
hand we have A\ = \/0.21. In general, for A e (x) — A(x/e) where A(y) = a + 



2G 



Figure 2: Example One: Here we have A\(x) plotted with the constant Ai(x) = 
V02T. 



/3 sin 2isy we have that A = a 2 — ft 2 . A plot of A\{x) is shown in Fig.jijand /3 is 
constant /3 = 0.01078280318. The parameters to this problem: T max = 12.4976 
(2000 HMM steps) and for all solvers except the exact finite difference solutions 
(direct numerical simulation, e.g., DNS) we used A = 0.5, p = 80. For DNS we 
used A = 0.5 and p e = 64. See Fig. [3] for a plot of the numerical solutions. 

4.2 Wave Propagation Problem Two 

We consider A £ (x) = 1.1 + \ (cos27ra + sin27ra;/e). The homogenized A can be 

computed analytically by freezing x: A{x) = ( A( ^ x ^ dyj . We have that 

A 2 (x) = v /(0.5cos(2 7 ra;) + l.l) 2 - 0.25. The profile A 2 (x) and /3 2 (x) is shown 
in Fig. U 

Numerical parameters: We have T max = 7.9985 (1600 time steps for HMM), 
A = 0.5 and p — 100 for all solvers except the DNS solver which uses A = 0.25 
and p e — 64. See Fig. [5] for a plot of the numerical solutions. 

4.3 Wave Propagation Problem Three 

Finally we consider A%{x) = (l.l + sin27r|) (1.5 + 0.5(cos(27ra; — 1))) which has 
the homogenized A 3 (x) = v / 0.21(1.5 + 0.5(cos(27ra;) - 1)). The profile A 3 (x) and 
f3a(x) is shown in Fig. [6} 

Numerical parameters: We have e = 0.03 and T max = 8.5677 (1200 HMM 
steps). For all sovers we used A = 0.5 and p = 70 except for the DNS solver 
which used A = 0.25 and p £ — 64. See Fig.[7|for a plot of the numerical solutions. 

In conclusion: We see that our numerical method accurately captures the 
solution of the effective equation. The choice of a bigger e is such that the 
model error in some cases is quite big. The difference between the exact and 
HMM solution and the solution of the effective equation will vanish as e — > 0. 
These illustrations show that our HMM method gives comparable results to the 
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Figure 4: Example Two: To the left we have A^x) plotted with A<z{x) and to 
the right we see (3i{x). 
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T=8.5714 




solution of an effective equation when it is known. But as we argued before; the 
effective equation might not be known. 



4.4 Numerical Investigation of the HMM Flux 

We will now demonstrate why the correction is essential to obtain a useful 
approximation to the long time wave propagation problem. We assume the 
standard wave propagation problem ([I]) where A\(x) — 1.1 + sin27r| and e 
0.01. We can then compute the correction matrix M and the uncorrected flux 
F, using the same micro solver parameters t = r\ — 20e, p E — 64, A = 0.5 and 
kernel from K 9,9 and five sample points. We get [\{M — /)!/] 



-2.4425e-ll 
-1.3116e-09 
5.2599e-07 
1.4890c-05 



-2.3810c-09 
1.5765e-10 
3.2966o-04 

-1.9744c-05 



4.7060c-02 
-3.3997e-08 
7.4096c-09 
2.0390c-03 



-2.8439c-10 
1.4118e-01 
3.9335e-05 

-2.3307c-06 



The corrected flux / is then (c.f. (25 1), 



2.303053047621222e-14 
4.582575696102362e-01 
9.160534264502929e-18 
1. 293944940509046o-05 



(44) 



/ = M f 



-2.083614204143469c-14 
4.582575696102412e-01 
2.386730479027304c-13 
6.469733186662828c-06 



(45) 



and the exact flux f e 



is 





V02l 


&pe 2 



p = 0.01078280318.. 



(46) 



o 

4.58257569495584e-01 


6.46968190800000c-06 



As we can see the correction affects mostly the last component of / related to 
the dispersive effects and not so much the lower order terms. The relative error 



f exact 



IS 



(fi - / x exact )//f xact is a mere 2.5019e-10 but the error (/ 3 - / 3 cxact )//; 
1.0000, i.e. zero correct digits. 

To conclude our investigation of the HMM fluxes we consider the second 
material A%{x) = 1.1 + 4 (cos 27ra; + sin ^p) which features both fast oscilla- 
tions and slowly varying parts. We will visualize the correction in the HMM 
procedure, 

\Jj__jj] Kj<3 



for the components 1, 2 and 3 corresponding to Q(x) — x l , 1 < i < 3 in the 
macroscopic points x = and x = 0.3. We also check the difference between 



the HMM corrected flux and the flux of the effective equation (43) where /3(x 
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was computed by freezing the slow variation in A\, i.e. 



\h-A{x)\ \U l/3-6 £ 2 /3(x)| 
e 2 e 2 ' e 2 

We let e vary between 0.1 and 0.001 in finite steps. We scale ij and r such that 
tj = t — 20e. In this experiment we used 16 points per e for the the micro solver, 
e.g., p £ = 16. See Figs.|8j[T0]and[T2]for plots of the convergence behavior. Note 
that /o and /o are both zero and therefore not plotted. 




Figure 8: A log-log plot of the difference 1 2 (left) and ^-5 — - (right) at 
x = for the material A\, using kernels with various p- and g-parameters. 




Figure 9: A log-log plot of the difference 1 2 1 (left) and -^-^ — - (right) at 
x = 0.3 for the material using kernels with various p- and g-parameters. 



We note that as e — > 0, compared to the effective flux there is an e 2 -difference 
in the components /1 and, when x = 0.3, also to /2- There is no difference in 
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Figure 10: A log-log plot of the difference ' ~ 2 ' (left) and '-j^ (right) at x = 
for the material Af , using kernels with various p- and g-parameters. Note that 
the effective equation flux is zero. 




Figure 11: A log-log plot of the difference lJ ^/^ (left) and ^ (right) at x = 0.3 
for the material A\, using kernels with various p- and g-parameters. Note that 
the effective equation flux is zero. 



fs and, when x = 0, no difference in fa. This gives us a clue to the form of 
the macroscopic flux in the case of a slowly variable coefficients A(x,x/e). We 
consider three possible candidates, 

^aiti = A.{x)u x + e 2 /3(x)u xxx , 
Fsiw = A(x)u x + e 2 d x j3(x)u xx , 
Faxts = A(x)u x + s 2 d xx /3(x)u x , 

where F a i t i is the effective equation used above. When u = x — x' where x' is 
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Figure 12: A log-log plot of the difference (left) and 1/3 ~jr 1 (right) at 

x = for the material Af, using kernels with various p- and g-parameters. 




the macroscopic point of evaluation we get 

F altl (x') = A(x'), 
F aU 2(x') = A(x'), 
F alt3 (x') = A(x')+e 2 f3"(x / ). 

Since we have a difference between /1 and F a \ t i(x') = A(x') of size e 2 it would 
indicate that the last alternative -Fkits is correct. Moreover, when u = (x — x') 2 , 
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we get 



F a iu(x') = 0, 
F alt2 {x') = 2e 2 f3'(x'), 
F a}t3 (x')=4e 2 /3'(x'). 

The e 2 -difference in J2 at x' = 0.3 but not at x' = is consistent with both 
F alt2 and F a i t3 , since /3'(0) = but /3'(0.3) ^ 0, cf. Fig g] Finally, when 
M = — a;') 3 then all three alternatives have the same value of 6e 2 f3(x') and 
the zero difference in / 3 is thus consistent with all of them. 

These results then lead us to tentatively conjecture the form of the flux F 
for the effective equation when the oscillatory coefficient has a slowly varying 
part: 

F = A(x)u x + e 2 d xx /3(x)u x . 

This would be the flux that our HMM procedure approximates. It differs from 
the effective equation ( 43 ) . Still the numerical results in Fig. p3]shows that HMM 



and the effective equation are very close in this case. The e -differences in the 
/1 and components seem to have very little effect on the solution over the 
time interval considered. On the other hand, the addition of the ^-component, 
which is also of size e 2 , has a dramatic effect. Not adding it corresponds to 



solving the homogenized equation ( 42 ) , which is far off the other solutions. We 
conjecture that for simulations over longer time the the differences in the /1 and 
ji components will play a bigger role. 



5 Conclusions 

We have developed and analyzed numerical methods for multiscale wave equa- 
tions with oscillatory coefficients. The methods are based on the framework of 
the heterogeneous multiscale method (HMM) and have substantially lower com- 
putational complexity than standard discretization algorithms. Convergence is 
proven in [9] for finite time approximation in the case of periodic coefficients and 
for multiple dimensions. The effective equation for long time is different from 
the finite time homogenized equation. After long time, dispersive effects enter 
and the method has to capture additional effects on the order of 0(e 2 ) [T5] . 
Numerical experiments show that the new techniques both accurately and ef- 
ficiently captures the macroscopic behavior for both finite and long time. It is 
emphasized that the HMM approach with just minor modifications accurately 
captures these dispersive phenomena. We prove that our method is stable if the 
spatial grid in the macro solver is sufficiently coarse. 
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